Mathematical modeling of circadian rhythms

Mathematical modeling of circadian rhythms#

circStudio implements four light-informed mathematical models designed to characterize and predict circadian physiology: Forger, Jewett, HannaySP and HannayTP. All four models can be used to estimate core body temperature minimum and the dim light melatonin onset. We begin by importing the necessary libraries:

import circStudio
import os

Import sample RPX file:

# Create file path for sample files within the circStudio package
fpath = os.path.join(os.path.dirname(circStudio.__file__))

# Create a new Raw instance using awd RPX adaptor
raw = circStudio.io.read_rpx(os.path.join(fpath, 'data', 'test_sample_rpx_eng.csv'), 
                             start_time='2015-07-04 12:00:00',
                             light_mode='White Light',
                             period='6 days',
                             language='ENG_UK')

1. Dynamic modeling of actigraphy-based light intensity data#

After loading the sample data, the next step is to create an instance of a circadian model for analysis. The model can be initialized using either a pd.Series or light data indexed by datetime, or two separate np.array objects, one for time and one for light intensity.

# Create a new instance of the HannaySP model using light data
model = circStudio.HannaySP(data=raw.light)

In this tutorial, we will demonstrate the functionality mainly using the HannaySP model, but syntax is consistent across the other available models:

# Jewett model
model = circStudio.Jewett(data=raw.light)

# Forger model
model = circStudio.Forger(data=raw.light)

# HannayTP model
model = circStudio.HannayTP(data=raw.light)

The resulting model states are stored in the model_states attribute. Be default, when a model is instantiated, the trajectory is computed using a predefined set of initial conditions, assuming a light schedule of 16 hours of light followed by 8 hours of complete darkness. To obtain more realistic and representative initial conditions for the model states, the model can be run iteratively over the light vector util it reaches a stable solution, assessed by convergence of the DLMO. For this, users must specify the number of iterations (loop_number) and provide the light data to be used in the loop:

# Loop over the provided light trajectory to determine new initial state
model.get_initial_conditions(data=raw.light, loop_number=50)

The model entrained after 6 loops.

array([  0.77016547, 263.37108225,   0.26764411])

From the model output, it is possible to extract a series of predicted DLMOs. By default, these values are expressed in hours relative to the start of the recording. To convert them to a 24-hour clock format, the modulo operator (% 24) can be applied:

# Obtain array of daily DLMOs
model.dlmos() % 24

array([7.525     , 7.4       , 7.28333333, 6.98333333, 6.7       ,
       7.3       ])

Using the plot() method in the Model class, it is possible to visually inspect the evolution of the daily predicted DLMO:

model.plot(dlmo=True)

Using the same function, it is possible to visualize all model states:

model.plot(states=True)

A key advantage of the HannaySP model, compared to the other available models, is its ability to compute the Entrainment Signal Regularity Index (ESRI). This metric quantifies how effectively a given light input can entrain the circadian system, ranging from 0 (no entrainment) to 1 (a maximally entraining signal).

# Compute ESRI
esri = circStudio.ESRI(data=raw.light)

# ESRI mean
esri.mean

np.float64(0.595468534541619)

# ESRI std
esri.std

np.float64(0.054790277969275325)

Future versions of circStudio will include additional mathematical models of circadian rhythms, as well as a model comparison utility. This functionality will allow users to systematically compare Kronauer-based, Van der Pol-type models with the physiologically-informed Hannay formalism. Stay tuned for upcoming developments!

2. Cosinor analysis#

# Create a cosinor object
cosinor = circStudio.Cosinor()

# Fit cosinor to raw activity data
results = cosinor.fit(data=raw.activity, verbose=False)

You can access individual best fit parameter values using the following syntax:

# Mesor
mesor = results.params['Mesor'].value

# Acrophase
acrophase = results.params['Acrophase'].value

# Period
period = results.params['Period'].value

# Amplitude
amplitude = results.params['Amplitude'].value

# Goodness-of-fit metrics (Akaike information criterium and Reduced Chi^2)
aic = results.aic
redchi = results.redchi

Additionally, you can also retrieve best fit parameters as a dictionary:

results.params.valuesdict()

{'Amplitude': np.float64(30.32292038752961),
 'Acrophase': np.float64(3.6510311570819654),
 'Period': np.float64(1626.5421520119878),
 'Mesor': np.float64(168.64347298661386)}

Finally, you can generate an interactive plot using the plot() method from Cosinor class:

cosinor.plot(data=raw.activity, best_params=results.params)

3. Next steps#

In the current version of circStudio (v1.0), support for Functional Linear Modeling (FLM), (Multifractal) Detrended Fluctuation Analysis (MF-DFA), and Singular Spectrum Analysis (SSA) has been temporarily removed to allow for reimplementation and enhanced functionality. These methods will be reintroduced in future versions of circStudio. Again, stay tuned!

Contents